str(Datafile_BP)
install.packages("metafor")
library(metafor)
install.packages("dplyr")
library(dplyr)

?escalc()

lnRR_BPa <- escalc(measure = "ROM", n1i = Datafile_BP$Samplesize_AGB_drought, n2i = Datafile_BP$Samplesize_AGB_control, m1i = Datafile_BP$AGB_drought, 
                  m2 = Datafile_BP$AGB_control, sd1i = Datafile_BP$SD_AGB_drought, sd2i = Datafile_BP$SD_AGB_control)
#where n1i and n2i are the sample sizes, m1i and m2i are the means, and sd1i and sd2i the standard deviations from each study
#The object created now has an effect size (yi) and its variance (vi) for each study
lnRR_BPt <- escalc(measure = "ROM", n1i = Datafile_BP$Samplesize_TB_drought, n2i = Datafile_BP$Samplesize_TB_control, m1i = Datafile_BP$TB_drought, 
                  m2 = Datafile_BP$TB_control, sd1i = Datafile_BP$SD_TB_drought, sd2i = Datafile_BP$SD_TB_control)
lnRR_BPb <- escalc(measure = "ROM", n1i = Datafile_BP$Samplesize_BGB_drought, n2i = Datafile_BP$Samplesize_BGB_control, m1i = Datafile_BP$BGB_drought, 
                  m2 = Datafile_BP$BGB_control, sd1i = Datafile_BP$SD_BGB_drought, sd2i = Datafile_BP$SD_BGB_control)


#now add the effect sizes to the original datafile
data_BPa <- bind_cols(Datafile_BP, lnRR_BPa)
data_BPt <- bind_cols(Datafile_BP, lnRR_BPt)
data_BPb <- bind_cols(Datafile_BP, lnRR_BPb)


#visualize the effect sizes with their conf intervals
forest(data_BPa$yi, data_BPa$vi)
#first 12 studies:
forest(data_BP$yi[1:12], data_BP$vi[1:12])

#random-effects model for agb

random_m_a <- rma(yi = yi, vi = vi, method = "REML", data = data_BPa)
#this gave a warning: Warning message: "Studies with NAs omitted from model fitting."

summary(random_m_a)


estimate = exp(-0.3372)          
#gives 0.7137661 = estimate (overall mean) of response ratios 
#statistically significant because p <0.0001
#reduction of 92% in aboveground biomass (0.71% when transformed back from ln: plants in drought treatments are 0.71% from the control biomass, so they are 29% smaller than the control)


#random-effects model for tb

random_m_t <- rma(yi = yi, vi = vi, method = "REML", data = data_BPt)
#this gave a warning: Warning message: "Studies with NAs omitted from model fitting."

summary(random_m_t)
estimate = exp(-0.7196)  

#gives 0.48695 = estimate (overall mean) of response ratios 
#statistically significant because p <0.0001

#random-effects model for bgb

random_m_b <- rma(yi = yi, vi = vi, method = "REML", data = data_BPb)
#this gave a warning: Warning message: "Studies with NAs omitted from model fitting."

summary(random_m_b)
estimate = exp(-0.3331)  

#gives 0.71667 = estimate (overall mean) of response ratios 
#statistically significant because p = 0.0068


#Meta-regression (stop predictors moderators in je model om heterogeniteit te verklaren)

time = data_BPa$'Drought_duration_(day)' 
duration = data_BPa$Duration
method = data_BPa$Experiment_method
fungrp = data_BPa$Group


metareg_a <- rma(yi = yi, vi = vi, mod = ~fungrp, method = "REML", 
               data = data_BPa)

#again the warning: "Studies with NAs omitted from model fitting"
summary(metareg_a)

inf <- influence(random_m_a)
print(inf)
plot(inf)

#R^2 = 0.00%; QM = 0.4055 --> they don't explain the heterogeneity

metareg_t <- rma(yi = yi, vi = vi, mod = ~fungrp, method = "REML", 
                 data = data_BPt)

summary(metareg_t)

# Egger's test (publication bias)
funnel(random_m)

# Note that the original Egger's test is regtest(random_m, model='lm')
regtest(random_m)
ranktest(random_m)

# Note that we are using the defult estimator ('L0'), but there are two
# others availablere
tf_m <- trimfill(random_m)
tf_m

funnel(tf_m)


#dplyr filter function --> use this to filter out the NAs and to only have true values
#tabel maken met alle aboveground biomass waarden zonder NA's


agb_noNA <- data_BPa%>%filter(AGB_drought > 0 & AGB_control > 0 & SD_AGB_drought > 0 & SD_AGB_control > 0
                                 & Samplesize_AGB_drought >0 & Samplesize_AGB_control >0)

agb_new <- select(agb_noNA, Paper_name, Paper_ID, yi, vi, Plant_species, Species_ID, Group, Duration, Experiment_type, AGB_drought, AGB_control, SD_AGB_drought, SD_AGB_control, Samplesize_AGB_drought, Samplesize_AGB_control, `Drought_duration_(day)`)
View(agb_new)

tb_noNA <- data_BPt%>%filter(TB_drought > 0 & TB_control > 0 & SD_TB_drought > 0 & SD_TB_control > 0
                             & Samplesize_TB_drought >0 & Samplesize_TB_control >0)

tb_new <- select(tb_noNA, Paper_name, Paper_ID, yi, vi, Plant_species, Species_ID, Group, Duration, TB_drought, TB_control, SD_TB_drought, SD_TB_control, Samplesize_TB_drought, Samplesize_TB_control, `Drought_duration_(day)`)
View(tb_new)

bgb_noNA <- data_BPb%>%filter(BGB_drought > 0 & BGB_control > 0 & SD_BGB_drought > 0 & SD_BGB_control > 0
                             & Samplesize_BGB_drought >0 & Samplesize_BGB_control >0)

bgb_new <- select(bgb_noNA, Paper_name, Paper_ID, yi, vi, Plant_species, Species_ID, Group, Duration, BGB_drought, BGB_control, SD_BGB_drought, SD_BGB_control, Samplesize_BGB_drought, Samplesize_BGB_control, `Drought_duration_(day)`)
View(bgb_new)

boxplot(formula = AGB_drought~`Drought_duration_(day)`, data = agb_new)

boxplot(formula = AGB_drought~Duration, data = agb_new, main="Effect of duration on biomass",
        xlab="Duration drought treatment", ylab="Above ground biomass (g)")

boxplot(formula = AGB_drought~Group, data = agb_new, main="Effect of functional group on biomass",
        xlab="Functional group", ylab="Above ground biomass (g)")

boxplot(formula = yi~Species_ID, data = agb_new, main="Effect sizes per species",
        xlab="Species", ylab="Effect size")

forest(agb_new$yi[1:25], agb_new$vi[1:25], xlab = "Effect size", main = "Studies aboveground biomass", header = "Distribution of effect sizes for aboveground biomass")

forest(agb_new$yi, agb_new$vi, xlab = "Effect size", header = "Distribution of effect sizes for aboveground biomass")


barplot(Estimate~Biomass, data = Effectsize_biomass, xlab = " ", ylab = " Effect size" , ylim = c(-0.9, 0.3), names.arg=c("Aboveground biomass","Belowground biomass","Total biomass"), col=rgb(0.2,0.4,0.6)) 

install.packages("ggplot2")
library(ggplot2)

# graph van verschillende biomassa op effectsize
ggplot(Effectsize_biomass) +
  geom_bar( aes(x=Biomass, y=Estimate), stat = "identity", fill="skyblue", alpha=0.8) + geom_text(aes(x = Biomass , y = Estimate, label=Number), vjust=-8.0, color="white", size=5) +
  geom_errorbar( aes(x=Biomass, ymin=Ci.lb, ymax=Ci.ub), width=0.2, colour="black", alpha=0.9, size=0.7) + labs(title="Effect of biomass on drought effects", x=" ", y = "Effectsize") + theme_minimal()

    
                                                                                                                
ggplot(Effectsize_functionalgroup) + geom_bar(aes(x=Functional_Group, y=Estimate), stat = "identity", fill="skyblue", alpha=0.8) + 
geom_errorbar( aes(x=Functional_Group, ymin=Ci.lb, ymax=Ci.ub), width=0.2, colour="black", alpha=0.9, size=0.7) + 
labs(title="Effect of drought per functional group", x=" ", y = "Effectsize") + theme_minimal()

forest(agb_new)


metareg_a <- rma(yi = yi, vi = vi, mod = ~duration, method = "REML", 
                 data = data_BPa)

summary(metareg_a)
res <- rma(yi = yi, vi = vi, method = "REML", data = agb_new)
#this gave a warning: Warning message: "Studies with NAs omitted from model fitting."

summary(res)
### set up forest plot (with 2x2 table counts added; the 'rows' argument is
### used to specify in which rows the outcomes will be plotted)
forest(res) 

FunGroup <- rma(yi = Estimate, vi = SE, data = Effectsize_functionalgroup)
forest(FunGroup)


ggplot(Effectsize_duration) + geom_bar(aes(x=Duration, y=Estimate), stat = "identity", fill="skyblue", alpha=0.8) + 
  geom_errorbar( aes(x=Duration, ymin=ci.lb, ymax=ci.ub), width=0.2, colour="black", alpha=0.9, size=0.7) + 
  labs(title="Effect of duration on drought effects", x=" ", y = "Effectsize") + theme_minimal()

y = metabin(data_BPa$yi, data_BPa$vi, data_BPa$Group, studlab= data_BPa$Paper_name)


forest(agb_new$yi[1:25], agb_new$vi[1:25])
forest(agb_new$yi, agb_new$vi)

forest(agb_new$yi[1:25], agb_new$vi[1:25], xlab = "Effect size", main = "Studies aboveground biomass", header = "Distribution of effect sizes for aboveground biomass")


forest(agb_new, xlim=c(-16, 6), at=log(c(0.05, 0.25, 1, 4)), atransf=exp,
       ilab=cbind(agb_new$Paper_name, agb_new$AGB_drought, agb_new$AGB_control, agb_new$yi, agb_new$vi),
       ilab.xpos=c(-9.5,-8,-6,-4.5), cex=0.75, ylim=c(-1, 27), rows=c(3:4,9:15,20:23),
       xlab="Risk Ratio", mlab="", psize=1, header="Author(s) and Year")

df <- data.frame(label, mean, lower, upper)

# reverses the factor level ordering for labels after coord_flip()
df$label <- factor(df$label, levels=rev(df$label))


dataforest <- ggplot(data=agb_new, aes(x= agb_new$Paper_name, y=agb_new$yi, ymin=agb_new$vi - agb_new$yi, ymax=agb_new$vi + agb_new$yi)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2) +  # add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("Study") + ylab("Estimate (95% CI)") +
  theme_bw()  # use a white background
print(dataforest)

p = ggplot(data=agb_new,
           aes(x =  agb_new$yi, y = agb_new$Paper_name, xmin=agb_new$yi - agb_new$vi, xmax=agb_new$yi + agb_new$vi, fill = agb_new$Group)) +
  geom_pointrange(aes(col=Group)) +  
  geom_hline(aes(fill=Group),yintercept =-0.5, linetype=2) +   xlab("Estimate (95% Confidence Interval)")+ ylab("Study") +
  geom_errorbar(aes(xmin=agb_new$yi - agb_new$vi, xmax=agb_new$yi + agb_new$vi,col=Group),width=0.5,cex=1) +
  #facet_wrap(~agb_new$Paper_name,strip.position="left",scales = "free_y") +
  theme(plot.title=element_text(size=16,face="bold"),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.x=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"))
   
print(p)

str(agb_new)

# even opnieuw de codes om resultaten te schrijven

# random effects model voor agb
random <- rma(yi = yi, vi = vi, method = "REML", data = agb_new)
#this gave a warning: Warning message: "Studies with NAs omitted from model fitting."

summary(random)

#Meta-regression (stop predictors moderators in je model om heterogeniteit te verklaren)

Group = agb_new$Group

metareg <- rma(yi = yi, vi = vi, mod = ~Group-1, method = "REML", 
                 data = agb_new)

#again the warning: "Studies with NAs omitted from model fitting"
summary(metareg)

# the real results

#the good forestplot
new_data<-agb_new%>%
  group_by(Paper_name, Group)%>%
  summarise(mean=mean(yi), sd=sd(yi))

apatheme=theme_bw() +
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Times'),
        legend.position='above',
        strip.text.y = element_text(size=10, color = 'black', face = 'bold'),
        strip.background = element_rect(colour = 'black', fill='powderblue',size=1, linetype='solid'))


p=ggplot(new_data, aes(y=Paper_name, x=mean))+
  #Add data points and color them black
  geom_point(aes(colour=Group))+
  #add the CI error bars
  geom_errorbarh(aes(xmin=mean-sd, xmax=mean+sd, colour=Group),height=.4)+
  #Specify the limits of the x-axis and relabel it to something more meaningful
  scale_x_continuous(limits=c(-1.25,1.25), name='Effect size')+
  #Give y-axis a meaningful label
  ylab('Reference')+
  #Add a vertical dashed line indicating an effect size of zero, for reference
  geom_vline(xintercept=0, color='black', linetype='dashed')+
  #Create sub-plots (i.e., facets) based on levels of setting
  #And allow them to have their own unique axes (so authors don't redundantly repeat)
  facet_grid(Group~., scales= 'free', space="free") +
  #apply the APA theme
  apatheme + theme(axis.title = element_text(size = 13))

print(p)

#Meta-regression (stop predictors moderators in je model om heterogeniteit te verklaren)

Group = agb_new$Group

metareg <- rma(yi = yi, vi = vi, mod = ~Group-1, method = "REML", 
               data = agb_new)

#again the warning: "Studies with NAs omitted from model fitting"
summary(metareg)

# plot the functional group + effect size
ggplot(Effectsize_functionalgroup) + geom_bar(aes(x=Group, y=Estimate, colour=Group, fill=Group), stat = "identity", alpha=0.8) + 
  geom_errorbar( aes(x=Group, ymin=Ci.lb, ymax=Ci.ub), width=0.15, colour="black", alpha=0.6, size=0.5) + 
  labs(title="Effect of drought on the aboveground biomass\nof functional groups", x=" ", y = "Effect size (% change compared to control)") +  theme_bw() +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank()) +  theme(axis.text.x = element_text(size = 15)) + theme(axis.title = element_text(size = 13)) + theme(plot.title = element_text(size = 14))     


# plot the functional group + effect size in percentage
ggplot(Effectsize_functionalgroup_) + geom_bar(aes(x=Group, y=Estimate, colour=Group, fill=Group), stat = "identity", alpha=0.8) + 
  geom_errorbar( aes(x=Group, ymin=Ci.lb, ymax=Ci.ub), width=0.2, colour="black", alpha=0.9, size=0.7) + 
  labs(title="Effect of drought on the functional groups", x=" ", y = "Effect size (% change compared to control)") + theme_minimal() 

# metaregression for drought duration
metareg_a <- rma(yi = yi, vi = vi, mod = ~duration-1, method = "REML", 
                 data = data_BPa)

summary(metareg_a)

# plot the drought duration + effect size
ggplot(Effectsize_duration) + geom_bar(aes(x=Duration, y=Estimate), stat = "identity", fill="mediumaquamarine", alpha=0.8) + 
  geom_errorbar( aes(x=Duration, ymin=ci.lb, ymax=ci.ub), width=0.15, colour="black", alpha=0.6, size=0.5) + 
  labs(title="Influence of drought duration on drought effects", x=" ", y = "Effect size (% change compared to control group)") + theme_bw() +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank()) +  theme(axis.text.x = element_text(size = 15)) + theme(axis.title = element_text(size = 13)) + theme(plot.title = element_text(size = 16))     

# plot the drought duration + effect size in percentage
ggplot(Effectsize_durationpercentage) + geom_bar(aes(x=Duration, y=Estimate), stat = "identity", fill="rosybrown1", alpha=0.8) + 
  geom_errorbar( aes(x=Duration, ymin=ci.lb, ymax=ci.ub), width=0.2, colour="black", alpha=0.6, size=0.5) + 
  labs(title="Influence of duration on drought effects", x=" ", y = "Effect size (% change compared to control group") + theme_minimal()

# plot the biomass + effect size
ggplot(Effectsize_biomass) +
  geom_bar( aes(x=Biomass, y=Estimate), stat = "identity", fill="mediumaquamarine", alpha=0.8) +
  geom_errorbar( aes(x=Biomass, ymin=Ci.lb, ymax=Ci.ub), width=0.15, colour="black", alpha=0.6, size=0.5) + labs(title="Effect of drought on biomass types", x=" ", y = "Effect size (% change compared to control group)") + 
  theme_bw() +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank()) +  theme(axis.text.x = element_text(size = 15)) + theme(axis.title = element_text(size = 13)) + theme(plot.title = element_text(size = 16)) + scale_x_discrete(labels=c("Aboveground\nbiomass", "Belowground\nbiomass", "Total\nbiomass"))

# plot the experiment type + effect size
ggplot(Effectsize_experimenttype) + geom_bar(aes(x=Experiment, y=Estimate), stat = "identity", fill="mediumaquamarine", alpha=0.8) + 
  geom_errorbar( aes(x=Experiment, ymin=Ci.lb, ymax=Ci.ub), width=0.15, colour="black", alpha=0.6, size=0.5) + 
  labs(title="Influence of experiment type on drought effects", x=" ", y = "Effect size (% change compared to control)") +  theme_bw() +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank()) + theme(axis.text.x = element_text(size = 15)) + theme(axis.title = element_text(size = 13)) + theme(plot.title = element_text(size = 16))     


# plot the biomass + effect size in percentage
ggplot(Biomass_effectsizepercentage) +
  geom_bar( aes(x=Biomass, y=Estimate), stat = "identity", fill="rosybrown1", alpha=0.8) + geom_text(aes(x = Biomass , y = Estimate, label=Number), vjust=-8.0, color="white", size=5) +
  geom_errorbar( aes(x=Biomass, ymin=Ci.lb, ymax=Ci.ub), width=0.2, colour="black", alpha=0.6, size=0.5) + labs(title="Effect of drought on biomass", x=" ", y = "Effect size (% change compared to control group)") + theme_minimal()

res.ml1 <- rma.mv(yi, vi, mod = ~Group-1, random = ~Group | Paper_name, data=agb_new)
summary(res.ml1, digits=3)
res.ml2 <- rma.mv(yi, vi, mod = ~Group-1, data=agb_new)
summary(res.ml2, digits=3)
anova(res.ml1, res.ml2)

metareg_a <- rma(yi = yi, vi = vi, mod = ~Group-1, method = "REML", 
                 data = agb_new)
summary(metareg_a)
anova(res.ml, metareg_a)

random_m_a <- rma(yi = yi, vi = vi, method = "REML", data = agb_new)

summary(random_m_a)

baujat(random)

res.ml <- rma.mv(yi, vi, random = ~ 1 | Paper_name/Group, data=agb_new)
print(res.ml, digits=4)

mean(agb_new$vi)

multilevel_a <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), 
                       method = "REML", data = agb_new)
summary(multilevel_a)

# random effects model
random_m_a <- rma(yi = yi, vi = vi, method = "REML", data = agb_new)
summary(random_m_a)

# Good random-effects models
multilevel_a <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), 
                       method = "REML", data = agb_new)
summary(multilevel_a)

multilevel_b <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), 
                       method = "REML", data = bgb_new)
summary(multilevel_b)

multilevel_t <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), 
                       method = "REML", data = tb_new)
summary(multilevel_t)

exp(-0.2468)

# Good meta-regression groups
multilevel_a_group <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), mod = ~Group-1,
                       method = "REML", data = agb_new)
summary(multilevel_a_group)

# Good meta-regression duration
multilevel_a_duration <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), mod = ~Duration-1,
                             method = "REML", data = agb_new)
summary(multilevel_a_duration)

# Good meta-regression experiment type
multilevel_a_exp <- rma.mv(yi = yi, V = vi, random = list(~1 | Paper_name, ~1 | Species_ID), mod = ~Experiment_type-1,
                                method = "REML", data = agb_new)
summary(multilevel_a_exp)


funnel(multilevel_a)

# Note that the original Egger's test is regtest(random_m, model='lm')
regtest(multilevel_a)

regtest(random_m_a)
